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developed  SLEDGE  computer  code .  This  simulation  was  done  £or 
seven  source  gamma  energies  from  0.5  to  7.0  MeV.  The  transport 
was  carried  to  a  range  of  25  mean  free  paths.  Correlated 
sampling  was  used  to  estimate  Compton  electron  currents  more 
accurately.  The  theoretical  basis  of  the  SLEDGE  code  is  dis¬ 
cussed.  Analytic  approximations  to  the  Monte  Carlo  results 
are  described  in  detail,  including  (1)  the  energy  deposition 
buildup  factor  (both  free  field  and  ground  influenced). ,  (2) 
time  dependence  of  energy  deposition  as  a  function  of  both  range 
from  the  burst  and  nearness  to  the  ground,  (3)  ratio  of  radial 
Compton  electron  current  to  energy  deposition  rate,  and  (4) 
theta  Compton  electron  current  due  to  the  presence  of  the  ground. 
Comparisons  are  made  with  two  previous  EMP  driver  packages,  with 
good  agreement  for  cases  where  the  previous  packages  are  be¬ 
lieved  reliable. 
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FOREWORD 


The  nuclear  electromagnetic  pulse  (EMP)  is  a  transient  broadband 
electromagnetic  field  capable  of  deunaging  or  upsetting  electronic  equip¬ 
ment  on  the  tactical  nuclear  battlefield.  Tto  predict  the  signature  of 
the  EMP  generated  by  a  nuclear  burst,  it  is  necessary  to  determine  the 
EMP  drivers — time  and  space  varying  ionization  and  Compton  electron 
currents  in  the  nuclear  radiation  field  around  the  burst.  Once  the 
ionization  and  Compton  currents  are  specified,  it  is  usually  possible  to 
solve  some  form  of  Maxwell's  equations  for  the  SIP  generated. 

This  report  summarizes  the  development  of  a  software  package  speci¬ 
fying  EMP  drivers  produced  by  prompt  gamma  rays  emitted  by  a  near¬ 
surface  burst  in  the  air.  the  package  was  designed  for  use  with  the 
HEMP  EMP  environment  prediction  computer  code  used  at  the  Harry  Diamond 
Laboratories.  A  separate  report,  A  Near-Surface  Burst  EMP  Driver 
Package  for  Neutron-Induced  Sources,  HDL-TR-1930  (September  1980), 
summarizes  the  development  of  a  similar  software  package  for  EMP  drivers 
due  to  neutrons  and  neutron -induced  secondary  gamma  rays. 


Part  I  of  this  report  describes  the  development  of  the  Monte  Carlo 
gamma  ray  transport  computer  code  and  data  base  to  be  fitted.  Part  II 
describes  the  curve-fitting  exercise  and  compares  these  results  with 
earlier  results. 
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PART  X.  BMP  DRIVERS  FROM  GAMMA  TRANSPORT  BY  MONTE  CARLO 


1.  REQUIREMENTS  FOR  NEW  EMP  DRIVER  AND  GAMMA  TRANSPORT  DATA 

Electromagnetic  pulse  (EMP)  drivers  consist  of  one  component  arising 
from  prompt  gamma  radiation  and  a  second  component  arising  from  neutron 
radiation.  The  NEMP  computer  code*'2  obtains  the  solution  of  Maxwell's 
equations  for  EMP  from  a  near-surface  burst  in  the  air.  The  EMP  drivers 
must  be  incorporated  into  the  NEMP  code  as  smooth  fits  to  quasi-analytic 
or  Monte  Carlo  predictions.  Earlier  EMP  driver  fits  were  used  in  the 
LEMP  computer  code3  for  EMP  from  surface  bursts.  The  LEMP  prompt  gamma 
ray  EMP  drivers  were  used  in  the  NEMP  code  at  first,  but  a  quick-look 
Monte  Carlo  study  by  this  author  confirmed  several  deficiencies  in  these 
LEMP  drivers: 

a.  Incorrect  range  dependence  for  the  gamma  ray  buildup  factor  at 
more  than  about  five  gamma  ray  attenuation  lengths  penetration,  becoming 
exponentially  worse  at  deeper  penetrations 

b.  Incorrect  time  dependence  at  deeper  penetrations  (LEMP  drivers' 
time  dependence  is  independent  of  range  and  decreases  with  time  too 
rapidly  at  deeper  penetrations.) 

c.  Lack  of  a  ground  effect  on  the  buildup  factor  (The  ground  actu¬ 
ally  inhibits  the  increase  of  the  buildup  factor  with  range  for  observ¬ 
ers  near  the  ground.) 

d.  Lack  of  late  vertical  Compton  currents  induced  by  the  effect  of 
the  ground  on  gamma  transport 

It  was  therefore  decided  to  develop  newer,  more  realistic  drivers 
for  the  NEMP  code  that  would  overcome  these  deficiencies  and  also  ac¬ 
count  for  the  effect  of  a  nonzero  burst  height.  The  new  drivers  would 
be  valid  to  a  penetration  depth  of  at  least  20  gamma  ray  attenuation 
lengths  and  would  encompass  more  discrete  gamma  ray  energies  than  the 
LEMP  drivers  (which  included  0.5-,  1.5-,  and  5.0-MeV  gamma  ray 
sources).  Hie  new  drivers  would  use  the  same  LEMP  philosophy  of  using 
well-chosen  analytic  functionals  to  describe  time  and  space  dependence 
of  the  data,  in  preference  to  using  piecewise  polynomial  fits. 

* H .  J.  Longle y,  C.  L.  Longmire ,  and  K.  S.  Smith,  Development  of  NEMP 
(U),  Mission  Research  Corp.,  Santa  Barbara,  CA,  HDL-CR-7 5-00 1-1  (April 
1975).  ( SECRET— RESTRICTED  DATA) 

2H.  J.  Longley  and  K.  S.  Smith,  Developments  in  NEMP  for  1977  (U), 
Mission  Research  Corp.,  Santa  Barbara,  CA,  HDL-CR-77-0022-1  ( January 
1978).  (SECRET — RESTRICTED  DATA) 

%.  J.  Longley  and  C.  L.  Longmire,  Development  and  Testing  of  LEMP  1, 
Los  Alamos  Scientific  Laboratory,  NM,  LA-4346  (April  1970). 


2.  DEVELOPMENT  OF  SLEDGE,  A  MONTE  CARLO  GAMMA  TRANSPORT  CODE 


To  first  test  the  validity  of  the  LEMP  drivers  and  then  obtain 
needed,  more  detailed  predictions  of  EMP  drivers,  it  was  decided  to 
write  a  special-purpose  Monte  Carlo  code  for  gamma  ray  transport,  spe¬ 
cially  tailored  for  efficiency  in  calculating  the  data  needed.  The  code 
was  written  in  FORTRAN  IV  and  run  on  the  Harry  Diamond  Laboratories 
( HDL)  computer  (an  IBM  System/370  Model  168).  flie  code  consists  of 
several  thousand  FORTRAN  statements  and  includes  53  program  and  subpro¬ 
gram  units.  It  was  named  SLEDGE,  after  the  basic  tool  for  realigning 
existing  structures;  "SLEDGE"  is  also  an  acronym  for  "surface  and  low- 
altitude  EMP  drivers  due  to  gamma  emission."  An  overview  of  the  SLEDGE 
code  is  presented  here  with  additional  detail  on  innovations  or  other 
important  aspects. 

The  SLEDGE  code  models  time- dependent  gamma  ray  transport  from  a 
point  source  in  air-over- ground  cylindrical  geometry.  The  collision 
kernel  includes  Compton  scattering,  photoabsorption,  and  pair  produc¬ 
tion.  Time  dependence  is  maintained  in  retarded  time  from  the  source 
point  (the  nuclear  burst).  For  this  study,  the  gamma  ray  source  was 
monochromatic  in  energy  and  a  delta  function  in  time. 

The  gamma  ray  flux  is  sampled  to  obtain  the  desired  EMP  driver 
response  functions — electron  current  vector  and  energy  deposition 
rate.  Energy  deposition  and  radial  and  theta  components  of  the  electron 
current  vector  are  scored  within  volume  detectors  by  both  track-length 
and  collision  estimators  of  the  gamma  ray  flux  and  are  scored  at  the 
surfaces  of  volume  detectors  by  a  boundary  crossing  estimator  of  the 
gamma  ray  flux.  Optionally,  collision  estimator  scoring  can  be  done  in 
spherical  polar  coordinate  radial  bins  to  obtain  one-dimensional  ( 1-d) 
results.  Thus,  removal  of  the  ground  allows  spherically  symmetric  free- 
field  results  to  be  accumulated  in  the  1-d  radial  bins.  Inclusion  of 
the  ground  allows  azimuthally  symmetric,  ground- influenced  results  to  be 
accumulated  in  the  2-d  cylindrical  geometry  volume  detectors. 

The  random  walk  permits  exponential  transform  biasing,  angular  and 
energy  biasing  of  the  source,  and  directional  biasing  of  gamma  rays 
emerging  from  Compton  collisions.  Splitting  can  be  done  to  prevent 
growth  of  particle  weight  due  to  biasing  at  collisions  or  at  trans¬ 
ports.  Splitting  and  roulette  can  be  used  in  connection  with  detector 
importances  to  improve  statistics  in  desired  regions  (zone  splitting). 
Splitting  and  roulette  can  be  used  also  to  reduce  the  variance  of  gamma 
ray  weights  in  specifiable  detectors  (weight  balancing).  A  special 
technique  may  be  used  to  greatly  reduce  the  variance  of  the  ratio  of 
electron  current  to  energy  deposition  rate,  based  on  the  (early  time) 
correlation  between  the  electron  currents  and  the  energy  deposition. 
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2.1  Random  Walk  and  Estimator  Scoring 

After  reading  problem  input,  setting  up  time  and  space  meshes, 
and  initializing  necessary  tables,  the  SLEDGE  code  executes  a  batch  loop 
and  then  scans  the  results  of  each  batch  to  develop  grand  totals  and 
batch  statistics.  Each  batch  consists  of  a  particle  (gamma  ray)  history 
loop  followed  by  saving  the  total  results  for  the  batch.  A  typical 
production  run  of  30,000  histories  might  consist  of  20  batches  of  1500 
histories  each.  Each  particle  history  consists  of  identifying  a  parti¬ 
cle  and  all  its  subparticles  created  by  splitting  and  following  them 
through  successive  random  walks  until  extinction.  The  random  walk 
consists  of  these  steps: 

( 1 )  Generating  the  particle  or  subparticle 

(2)  Determining  the  number  of  mean  free  paths  down  which  the 
particle  is  to  be  transported  prior  to  a  collision 

(3)  Tracking  the  particle  through  space  and  time  to  the  colli¬ 
sion  point,  while  scoring  track- length  or  surface  crossing  estimators 
for  the  detectors  encountered 

(4)  Colliding  the  particle  to  obtain  the  emergent  particle 
energy  and  direction,  while  scoring  collision  estimators  for  the  detec¬ 
tor  involved,  and  then  repeating  steps  (2),  (3),  and  (4)  until  extinc¬ 
tion 


The  particle  or  subparticle  can  be  killed  or  split  after  step 
(2)  or  (3)  or  (4)  through  photoabsorption,  escape,  age  cutoff,  low 
energy  cutoff,  collision  number  limit,  or  any  of  the  roulette  or  split¬ 
ting  mechanisms.  Subparticles  produced  are  stored  in  a  split  stack  and 
are  withdrawn  one  at  a  time  from  the  split  stack  in  step  ( 1 ) .  Ultimate¬ 
ly,  there  are  no  entries  remaining  in  the  split  stack,  and  the  next 
particle  history  can  be  started. 

Many  standard  techniques  were  drawn  from  the  literature,  par¬ 
ticularly  from  Carter  and  Cashwell.1* 

Generating  the  particle  involves  determining  its  position  (at 
the  burst),  direction  cosines  (usually  biased),  energy  (simply  monochro¬ 
matic  in  this  study),  age,  and  weight. 

Determining  the  number  of  mean  free  paths,  N,  was  done  by 

N  =  -  in  (K) 


*♦ L .  L.  Carter  and  E.  D.  Cashwell,  Particle-Transport  Simulation  with  ' 
the  Monte  Carlo  Method ,  U.S.  Energy  Research  and  Development 

Administration  TID-26607  (1975) . 
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for  a  uniform  random  number,  £,  selected  from  the  interval  (0,1).  A 
rejection  method  was  used  to  eliminate  escapes  before  the  first  scatter 
for  a  particle.  Both  exponential  transform  and  importance  function 
biasing  were  used  for  first-scatter  transport,  and  exponential  transform 
biasing  was  used  for  later  scatters.  These  are  discussed  in  more  detail 
in  section  2.2. 

Tracking  the  particle  involved  finding  its  track-length  in  the 
detector  that  it  was  inside  (detector  volumes  being  mutually  exclusive 
and  filling  all  space  within  the  grid  outer  boundary) ,  scoring  track- 
length  and  surface  crossing  estimators,  moving  the  particle  to  the  next 
detector,  and  repeating  the  process  until  the  particle  reached  a  col¬ 
lision  point  or  was  killed.  When  the  particle  reached  a  collision 
point,  the  2-d  and  1-d  collision  estimators  were  scored,  also. 

Standard  means  were  used  to  select  the  collision  event  from 
Compton  collision,  photoabsorption,  and  pair  production  probabilities. 
Pair  production  was  assumed  isotropic.  Standard  algorithms  were  used 
to  select  the  emergent  energy  and  the  direction  of  Compton  scattered 
gamma  rays.  Emergent  direction  biasing  could  be  used,  as  described  in 
section  2.2. 

The  track-length  estimation  of  the  flux,  F  (that  is,  path 
length  per  unit  volume),  is 


F  -  pWBC/V  , 

where  p  is  the  track-length  in  the  particular  volume  detector  and  time 
bin,  Woc  is  the  particle  weight  before  collision,  and  v  is  the  detector 
volume.  The  electron  current  magnitude,  J,  is 

J  =  F  (Kc/Xc  +  RpAp)  ' 

where  Rc  is  the  average  forward  range  of  Compton  electrons  ( for  the 
incident  gamma  ray  energy) ,  Xc  is  the  Compton  collision  or  attenuation 
length,  Rp  is  the  average  forward  range  of  photoelectrons,  and  Xp  is 
the  photoelectric  absorption  length.  Rc  was  computed  and  tabulated  by 
the  SI£DGE  code  by  using  the  Compton  differential  scattering  cross 
section  and  the  electron  range-energy  relation  in  air.  Rp  was  adapted 
from  the  literature.  The  energy  deposition,  D,  is 

D  =  F[B/Xa  +  E/Xp  +  (E  -  1 • 022 )/XppJ  , 

where  E  is  the  gamma  ray  energy  in  MeV,  Xg  is  the  Compton  absorption 
length,  and  Xpp  is  the  pair  production  absorption  length. 

Once  J  is  obtained,  it  is  decomposed  into  radial  and  vertical 
components • 
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The  surface  crossing  estimation  of  the  flux,  scored  only  when 
the  surface  is  crossed,  is 

F  =  Wbc/CA  , 

where  C  is  the  sine  of  the  angle  between  the  surface  and  the  particle 
direction  and  A  is  the  area  of  the  surface.  J  and  D  are  derived  as 
before. 

The  collision  estimation  of  the  flux,  scored  only  when  a  col¬ 
lision  point  is  reached,  is 

f  =  wBC/<rv  , 

where  a  is  the  total  collision  cross  section  in  the  detector  with  volume 
V.  J  and  D  are  derived  as  before. 

The  time  bin  associated  with  scoring  surface  crossing  and 
collision  estimators  is  simple  to  determine.  However,  for  the  track- 
length  estimator,  the  particle  may  pass  through  several  time  bins  while 
traversing  a  detector  volume.  The  portion  of  the  track- length  spent  in 
each  time  bin  must  be  determined  to  score  the  time  dependence  prop¬ 
erly.  Although  this  determination  requires  somewhat  more  computer  time, 
the  associated  variance  is  usually  much  better  for  the  track-length 
estimator.  Generally,  all  three  estimators  were  used  to  check  for 
agreement  and  insure  correct  results. 

No  account  was  taken  of  the  effect  of  the  ground  upon  the 
transport  of  Compton  electrons  and  photoelectrons.  Some  effects  would 
be  expected  within  roughly  1  m  of  the  ground.  This  effect  could  be  a 
fruitful  subject  for  future  study. 

2.2  Biasing  and  Importance  Functions 

The  analog  model  of  gamma  ray  transport  (that  is,  no  biasing  or 
importance  functions)  is  useless  for  deep  penetration  problems.  Some 
means  must  be  used  to  force  gamma  rays  to  penetrate  to  distant  regions 
in  large  enough  numbers  to  produce  useful,  low  variance  results.  A 
number  of  methods  were  tested  in  the  SLEDGE  code  for  the  problems  under 
consideration.  These  may  be  categorized  as  (1)  source  biasing,  (2) 
transport  biasing,  and  (3)  collision  biasing,  modifying  respectively  the 
source  function,  the  transport  kernel,  or  the  collision  kernel  of  the 
transport  integral.  Two  risks  are  associated  with  the  use  of  properly 
normalized  biasing  methods: 

a.  The  variance  of  the  flux  may  be  infinite  or  possibly  very 
large  so  that  as  the  number  of  batches  increases,  the  variance  gets 
larger,  not  smaller. 


b.  The  transport  integral  (integrated  over  position,  velocity, 
time,  and  energy)  may  be  undersampled  since  strong  biasing  may  force  the 
sampling  away  from  portions  of  the  integral  space  that  contribute  sig¬ 
nificantly  to  the  correct  answer. 

Thus,  the  biasing  methods  must  be  properly  constructed  and  not 
overused  and  yet  must  be  powerful  enough  to  be  worth  the  trouble  of 
using.  In  general,  source  biasing  very  effectively  reduces  the  variance 
in  deep  penetration  problems.  Also,  splitting  with  roulette  is  a  simple 
but  highly .  effective  method  for  transport  biasing.  In  this  study, 
collision  biasing  was  only  of  minor  importance  in  variance  reduction. 

2.2.1  Source  Biasing 

Source  directional  biasing  was  not  used  in  the  SLEDGE  code  to 
obtain  1-d  results.  However,  for  the  2-d  air-over-ground  studies, 
source  directional  biasing  was  used  by  construction  of  a  cumulative 
distribution  function  (cdf)  based  on  a  nonisotropic  importance  func¬ 
tion.  Conventional  cylindrical  coordinates  (p,z,$)  are  used,  where  z  is 
positive  upward.  If  the  particle  direction  cosine,  x,  is 

x  =  z/(p2  +  z2)l/2  , 

then  the  probability  density  function  (pdf) ,  which  describes  an  iso¬ 
tropic  source  angular  distribution,  is 

p(x)  =  1/2,  -1  <  x  <  1  . 

A  suitable  anisotropic  pdf,  which  biases  the  selection  toward  certain 
target  direction  cosines  x^,  is 

N  , 

P ( x)  =  A  l  jl  +  Bj/rl  +  B^x  -  xi 
i=1  X 


where  A  normalizes  p(x)  to  have  an  integral  of  unity  (taken  from  x  =  -1 
to  x  **  1)  and  where  3^  and  x^  are  N  bias  strengths  and  bias  direc¬ 
tions.  If  a  is  zero,  there  is  no  bias  toward  the  corresponding  x^ ; 
if  all  6^  are  zero,  the  pdf  is  isotropic.  The  associated  cdf  is 


c(x) 


dx"  . 


and  the  scattering  direction  is  obtained  in  the  usual  way  by  selecting  a 
random  number,  £ ,  from  the  uniform  distribution  on  (0,1)  and  then  sol¬ 
ving 

c(x)  =  £ 


for  x.  The  associated  weight  is  simply  1/p(x) 


The  transport  to  the  £irst  collision  was  included  as  part  of 
the  source  biasing  and  treated  separately  from  subsequent  transport 
biasing.  In  addition  to  the  directional  biasing  already  described,  two 
types  of  path  length  biasing  were  used:  (1)  the  well-known  exponential 
transform  and  (2)  an  importance  function.  The  exponential  transform  is 
not  described  further  since  it  is  a  standard  technique.  For  the  impor¬ 
tance  function,  instead  of  sampling  from  a  simple  exponential  distribu¬ 
tion,  e~ar,  in  the  transport  integral 

/Se ~ar  dp  , 

where  a  is  the  range  attenuation  coefficient,  we  wish  to  sample  from  the 
distribution 


(1  +  Yr)e‘nr  , 

where  y  is  an  additional  parameter  and  n  is  a  different  attenuation 
coefficient.  Typically,  n  may  be  an  order  of  magnitude  smaller  than  a 
to  achieve  deep  penetration  by  the  first  scatter,  and  y  may  be  about 
unity.  If  y  is  zero,  the  importance  function  reduces  to  the  exponential 
transform,  while  if  y  is  about  unity,  the  deep  penetration  statistics 
are  markedly  improved  over  the  simple  exponential  transform.  Since 
gamma  ray  buildup  factors  are  roughly  proportional  to  r,  the  yr  factor 
actually  suffices  to  compensate  for  1/r2  geometric  attenuation.  A 
rejection  method  also  was  used  to  eliminate  escapes  before  the  first 
collision. 


2.2.2  Transport  Biasing 

Several  methods  of  biasing  the  transport  kernel  were  used.  A 
directional  exponential  transform  was  sometimes  used  to  force  longer 
transports  in  the  outer  direction,  but  was  found  unnecessary  when  region 
weight  balancing  and  region  importances  were  used  with  splitting  and 
roulette . 


Weight  balancing  involved  keeping  the  particle  weights  near 
an  optimum  weight,  which  was  a  function  of  the  distance  from  the 
burst.  As  necessary,  particle  weights  were  adjusted  after  each  col¬ 
lision  by  splitting  or  playing  roulette.  Keeping  the  weights  nearly  the 
same  within  a  volume  detector  caused  the  mean  square  weight  to  be  close 
to  the  mean  weight,  leading  to  improved  statistics.  The  optimum  weight 
was  the  first  collision  weight  multiplied  by  the  approximate  buildup 
factor, 

wopt  =  B(r)wKr)  , 

where  B(r)  is  the  approximate  buildup  factor  at  range  r  from  the  burst 
and  wi(r)  is  the  first-collision  weight  at  range  r.  If  no  source  trans¬ 
port  bias  were  used,  weight  wi(r)  of  a  particle  undergoing  a  first 
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collision  at  range  r  would  be  unity  for  all  r.  However,  the  use  of  the 
importance  function  described  in  section  2.2.1  causes  weight  wi  to  be  a 
function  of  r. 


Region  importances  also  were  implemented.  Various  regions 
were  assigned  numerical  ratings  of  their  importance  by  the  input  to  the 
SLEDGE  code.  For  example,  detectors  near  the  ground  were  more  important 
than  detectors  far  above  the  ground.  When  a  particle  passed  from  a 
region  of  less  importance  to  a  region  of  greater  importance,  it  was 
split.  When  passing  the  other  way,  roulette  was  played.  ffius,  more 
particles  were  tracked  in  more  important  regions,  and  fewer  were  tracked 
in  less  important  regions. 

2.2.3  Collision  Biasing 

In  several  problems,  it  was  found  helpful  to  use  collision 
biasing,  in  which  the  emergent  particle  direction  (scattering  direction) 
was  biased  toward  one  or  more  areas  of  particular  interest.  Die  scat¬ 
tering  direction  was  biased  only  for  Compton  scattering^  For  photo¬ 
absorption,  the  gamma  ray  is  completely  absorbed.  For  pair  production, 
the  scattering  directions  of  the  two  locally  produced  annihilation 
quanta  were  assumed  to  be  isotropic  and  opposite  to  each  other .  For 
Compton  scattering,  the  bias  importance  function  used  was 

*c  *  l  V[l-  <i(***i)]  * 

i=1  \  /J 

where 

is  the  bias  strength  in  the  ith  bias  direction, 

A 

r  is  the  unit  vector  in  the  scattering  direction, 

A 

r ^  is  the  unit  vector  in  the  ith  bias  direction. 

Hie  dot  product  r «r^  is  the  cosine  of  the  angle  between  r  and  r^  so  that 
Ic  is  maximized  when  r  and  r*  are  parallel. 

Construction  of  a  properly  normalized  pdf  based  on  importance 
function  Xc  cm  be  awkward,  because  the  complete  angular  integral, 

Q  mJ4ll  *  Vc«M)  ' 

is  required,  where  oc($, 0)  is  the  Compton  differential  scattering  cross 
section.  Since  oc($,6)  is  nkt  a  simple  function,  the  integral  cannot  be 
obtained  analytically.  A  numerical  integration  would  be  time- 
consuming.  Even  a  straightforward  rejection  method  requires  the  same 


integral  to  conserve  particle  weight.  Hence,  a  special  algorithm  was 
devised  to  allow  use  of  the  rejection  method  to  implement  Ic.  The 
algorithm  uses  recursive  binary  splitting  with  roulette  to  obtain  ex¬ 
actly  one  surviving  particle.  The  algorithm  conserves  weight  and  has  a 
low  variance.  Briefly,  the  algorithm  consists  of  the  following  steps  to 
obtain  the  biased  scattering  direction,  Q,  for  a  particle  of  weight,  W: 

(1)  Select  roulette  variable  V  uniformly  from  the  set  (1,2). 

(2)  Generate  random  trial  scattering  direction  S)i  for  a 

first  split  and  get  rejection  probability  Ri . 

(3)  Test  for  rejection  of  the  first  split.  If  it  is  re¬ 
jected,  set  W  =  W/2  and  return  to  step  (2). 

(4)  Otherwise,  generate  random  trial  scattering  direction 
fl2  for  a  second  split  and  get  rejection  probability  R2. 

(5)  Test  for  rejection  of  the  second  split.  If  it  is  re¬ 
jected,  set  W  =  W/( 2  -  2/Rl ) ,  set  (2  =  fti ,  and  go  to  Finish. 

(6)  Otherwise,  if  V  =*  1,  set  W  =  W/(1  -  Ri ) ,  set  ft  »  ftl, 

and  go  to  Finish. 

(7)  Otherwise,  if  V  =  2,  set  W  *  W/( 1  -  R2),  set  ft  *  02, 

and  go  to  Finish. 

Finish  This  is  the  end  of  the  algorithm. 

Ibis  algorithm  allows  the  use  of  the  rejection  method  without  the  need 
to  evaluate  normalization  integral  Q. 

In  a  number  of  SLEDGE  runs,  this  algorithm  was  successfully 
used  to  push  particles  away  from  the  ground  and  toward  certain  volume 
detectors  with  poor  statistics. 

2.3  Cross  Sections 


The  well-known  Klein-Nishina  expressions  were  used  for  the 
Compton  scattering  cross  section.  For  the  media  under  consideration 
(air  and  ground) ,  the  Compton  cross  section  dominates  cross  sections  for 
photoabsorption  and  pair  production  over  a  wide  energy  range. 


The  attenuation  coefficient  for  pair  production  was  approxi¬ 
mated  by  the  simple  expression 


PP 


1.84  x  10_4pZ2(e  -  0.014E2)/A, 


a 


where 


T _  is  in  inverse  centimeters, 

PP 

p  is  the  material  density  in  grams  per  cubic  centimeter, 

22  is  the  mean  square  atomic  number  of  the  material, 

a 

E  is  E  -  1.022,  and  E  is  the  gamma  ray  energy  in  MeV, 

A  is  the  mean  atomic  weight  of  the  material. 

For  air,  Z2  =  106.4  and  A  =  28.9. 

For  aluminum,  Z2  =  132  and  A  =  27.0. 

The  attenuation  coefficient  for  photoabsorption  was  approxi¬ 
mated  by 

1.88  x  10“2pZ4 

Tph  ’  A<100E)3-49-0*0035(1o9  54)2  ' 


where  the  logarithm  is  to  the  base  10.  1116  material  mean  fourth  power 

of  the  atomic  number  is  6.46  x  lo3  for  air  and  131*  for  aluminum.  This 
approximation  has  no  significant  error  for  air,  but  is  about  10  percent 
low  for  aluminum. 


The  error  incurred  in  the  total  absorption  cross  section  by 
using  these  approximations  appears  to  be  less  than  1  percent  in  air  for 
gamma  ray  energies  from  0.01  to  10  MeV.  The  error  in  aluminum  is 
slightly  greater  for  energies  less  than  0.1  MeV.  The  values  for  alumi¬ 
num  are  cited  here  because  the  air-over- ground  transport  studies  were 
done  over  an  aluminum  ground  for  simplicity.  The  gamma  ray  cross- 
section  properties  of  aluminum  are  similar  to  a  typical  ground  composed 
of  silicon,  aluminum,  oxygen,  and  water. 

2.4  Scoring  EMP  Drivers  by  Correlated  Sampling 

A  particularly  useful  innovation  in  the  SLEDGE  code  concerns 
exploitation  of  the  correlation  between  the  energy  deposition  response 
function  and  the  (Compton  and  photoelectric)  electron  current  response 
function.  When  EMP  drivers  are  scored,  the  electron  current  contribu¬ 
tion  is  found  to  be  correlated  with  the  associated  energy  deposition 
contribution.  The  correlation  is  high  at  early  times  when  the  gamma  ray 
direction  is  well  collimated,  but  becomes  low  at  later  times  when  the 
gamma  ray  direction  is  poorly  collimated.  Ultimately,  the  correlation 
drops  to  zero  as  the  marfy  scatters  produce  am  isotropic  gamma  direc¬ 
tion.  The  correlation  can  be  exploited  in  the  following  way. 


Define  the  energy  deposition  score  in  a  volume  detector,  k,  and 
time  bin,  l, 


n 

skA  =  n  l  akto 
rn=1 

and  the  radial  ( for  example)  electron  current  as 

n 

Jkf  “  n  l  * 

m=1 

where  m  denotes  the  batch  number  of  a  total  of  n  batches.  Thus,  is 

the  energy  deposition  score  for  the  kth  detector,  ith  time  bin,  and  mth 
batch.  Define  also  arbitrary  function  which  depends  on  the  detec¬ 

tor  and  the  time  bin.  Consider  the  variance  of  the  linear  combination 

Jk*  ~  "kf^f5 

VAR  {Jkf  “  Hkfskf}  =  £  jkfm  ‘  2Hkf  l  jktoskto  +  a2kJl  l  skto  ‘ 

'  )  m  m  m 


If  and  are  positively  correlated,  the  negative  term  on  the 

right  hand  side  reduces  the  variance*  if  the  correlation  coefficient  of 
jkfm  and  Skim  is  unity,  the  right  hand  side  becomes  zero.  Therefore, 
on  the  assumption  that  jj^  and  are  at  least  partly  correlated, 

solve  the  equation 

jk*m  =  HkJlSkJlm 

for  by  least-squares  regression.  This  regression  value  for 

should  minimize 


VAR  [\t  " 

and  make  optimum  use  of  any  correlation  between  and  Hk£sk£m* 

In  practice  with  the  SLEDGE  code,  this  correlated  sampling 
technique  improved  the  electron  radial  current  fractional  deviation 
statistics  by  more  than  an  order  of  magnitude ,  except  at  later  times 
when  all  remaining  gamma  rays  have  been  scattered  many  times.  The 
electron  theta  current  statistics  were  similarly  improved  by  more  than  a 
factor  of  two,  especially  at  late  times,  when  significant  theta  currents 
were  developing. 
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3.  CALCULATIONS  PERFORMED 


3.1  One-Dimensional  Studies 


The  principal  aspects  of  gamma  ray  transport  can  be  described 
in  terms  of  buildup  factors.  (A  buildup  factor  for  an  effect  is  the 
ratio  of  the  effect  as  engendered  by  uncollided  and  multiply  scattered 
gamma  ray  flux,  to  the  effect  as  engendered  by  uncollided  gamma  flux.) 
Hie  SLEDGE  code  was  run  for  seven  discrete  gamma  ray  energy  sources 
(table  1)  in  a  homogeneous  free  air  geometry.  Results  for  the  time 
integrated  energy  deposition  buildup  factor  were  obtained  from  25  1-d 
radial  bins,  each  1  attenuation  length  wide,  extending  to  a  maximum 
range  of  25  attenuation  lengths.  Fractional  deviations  v.aried  from 
about  3  percent  at  the  closer  radial  bins  to  about  10  percent  for  the 
most  distant.  For  each  gamma  ray  energy,  a  SLEDGE  run  was  made  con¬ 
sisting  of  25  batches  of  4000  histories  each,  for  a  total  of  100,000 
histories.  A  typical  run  tallied  about  one  million  collisions  (scat¬ 
ters)  and  consumed  about  40  min  of  central  processing  unit  (CPU)  time  on 
the  HDL  computer  (an  IBM  370/168).  Results  were  stored  on  magnetic  tape 
for  subsequent  analysis  and  curve  fitting. 


TABLE  1.  SEVEN  GAMMA  RAY  ENERGIES  AND  ATTENUATION  PARAMETERS 


Source  gamma  ray 
energy 
(MeV) 

Total  attenuation 
coefficient 
(o_1 ) 

1  attenuation 
length 
(m) 

25  attenuation 
lengths 
(m) 

0.5 

0.009604 

104. 12 

2603.08 

1.0 

0.007014 

142.57 

3564.30 

1.5 

0.005734 

174.40 

4359.96 

2.5 

0.004377 

228.47 

5711.67 

3.5 

0.003651 

273.90 

6847.44 

5.0 

0.003034 

329.60 

8239.95 

7.0 

0.002596 

385.21 

9630.20 

Table  1  lists  the  total  attenuation  coefficient  and  distances 
corresponding  to  1  and  25  attenuation  lengths.  An  air  density  of 
1.11  x  kg/m3  is  assumed. 

Time-dependent  results  also  were  obtained  from  these  studies. 
The  time  bin  structure  used  was  the  same  as  for  the  2-d  problems,  shown 
in  table  2. 
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TABLE  2.  SLEDGE  MESHES  FOR  TWO- 
DIMENSIONAL  PROBLEMS 


Time  bin 
boundaries 
(s) 

Radial 

boundaries 

(m) 

Vertical 

boundaries 

(m) 

0.00 

0.0 

-1.0 

5.00(-9)a 

300.0 

•Q 

o 

o 

1.001-8) 

600.0 

25.0 

2.00 ( -8) 

900.0 

50.0 

3.251-8) 

1200.0 

100.0 

5.00 ( -8) 

1500.0 

200.0 

7. 50 (-8) 

1800.0 

400.0 

1 . 00 ( —7 ) 

2100.0 

600.0 

1 . 50 ( -7 ) 

2400.0 

3000.0 

2. 00( -7 ) 

2700.0 

3.00( -7) 

3000.0 

4.001-7) 

6.00(-7) 

1.00<-6) 

1.50 (-6) 

2.00<-6) 

3.00 ( -6) 

5.00<-6) 

1.001-5) 

aRead  as  5.00  x  I0”9. 
b Air-ground  interface  at  O.O- m  height. 

3.2  Two-Dimensional  Air-Ground  Interface  Effects 

The  2-d  air-over-ground  gamma  ray  transport  problem  may  be 
considered  as  a  perturbed  free-field  problem.  Significant  perturbations 
are  expected  near  the  air-ground  interface.  The  SLEDGE  code  was  run  for 
a  1.5-MeV  gamma  ray  source  in  a  cylindrical  air- over- ground  geometry.  A 
burst  height  of  200  m  was  used.  As  mentioned  in  section  2.3,  the  ground 
was  assumed  for  simplicity  to  be  alumimm,  whose  gamma  ray  cross-section 


properties  are  analogous  to  a  typical  ground  composed  of  silicon,  alumi¬ 
num,  oxygen,  and  water.  The  radial  (horizontal)  mesh  extended  to  a 
3000-m  ground  range,  and  the  vertical  mesh  extended  to  a  3000-m  height 
aboveground  and  to  a  1-m  depth  belowground.  Eighteen  time  bins  were 
used,  extending  to  a  10-ys  retarded  time.  lhese  grids  are  shown  in 
table  2.  Correlated  sampling  was  used  for  the  radial  Compton  current, 
but  not  for  the  theta  Compton  current.  Various  biases  and  importance 
functions  were  used,  as  described  in  section  2.2.  The  run  consisted  of 
30  batches  of  20,000  histories  each,  for  a  total  of  600,000  histories. 
In  longer  than  2  hr  of  CPU  running,  nearly  four  million  collisions 
(scatters)  were  tallied.  Fractional  deviations  for  energy  deposition 
averaged  1  percent  close  to  the  burst  and  about  6  percent  near  the  3-km 
boundary . 

It  was  observed  that  energy  deposition  from  outgoing  gamma  rays 
was  about  eight  times  that  from  backscattered  incoming  gamma  rays,  at 
any  range  in  the  air.  'Hie  backscatter  from  the  ground  was  much  less,  so 
that  radial  Compton  current  and  energy  deposition  were  weaker  near  the 
ground  and  the  theta  Compton  current  was  strongest.  This  effect  is 
clear  since  gamma  rays  must  flow  from  the  air  into  the  ground  and  be 
attenuated  there;  this  flow  causes  a  negative  gradient  in  the  flux  close 
to  the  ground.  From  a  different  viewpoint,  the  buildup  factor  is  re¬ 
duced  near  the  ground. 

Since  effects  due  to  multiply  scattered  gamma  rays  come  later 
in  time  than  effects  due  to  unscattered  gamma  rays,  the  time  dependence 
of  effects  was  confirmed  to  be  range  dependent.  At  deeper  penetrations, 
the  buildup  factor  is  greater  so  that  multiply  scattered  effects  are 
stronger.  Ttoe  result  is  that  the  energy  deposition  decay  rate  in  time 
is  slower  at  deeper  penetrations.  The  ratio  of  radial  Compton  current 
to  energy  deposition  rate  also  declines  with  time  more  slowly  with 
increasing  range.  In  fact,  the  time  dependence  of  the  latter  ratio 
roughly  scales  as  the  inverse  of  the  buildup  factor  (that  is,  scaled 
time  t'  «  tB) . 

A  similar  SLEDGE  run  was  done  for  a  0.5-MeV  gamma  ray  source. 
Results  of  both  runs  were  stored  on  magnetic  tape  for  subsequent  analy¬ 
sis. 
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PART  II.  ANALYTICAL  FITS  TO  EMP  DRIVER  RESULTS 


4.  ENERGY  DEPOSITION 

4.1  Fitting  Free-Field  Buildup  Factors  for  Seven  Gamma  Ray  Energies 

Seven  gamma  ray  energies  from  0.5  to  7.0  MeV  (table  1)  were 
chosen  to  allow  flexibility  in  representing  approximately  an  arbitrary 
nuclear  burst  gamma  ray  energy  spectrum.  Since  gamma-induced  EMP 
sources  may  contribute  to  EMP  generation  at  distances  as  great  as  about 
5  km,  the  gamma  ray  transport  was  carried  to  a  range  of  25  attenuation 
lengths  to  ascertain  the  true  range  dependence  of  the  buildup  factor  for 
energy  deposition.  Based  on  theory,  a  reasonable  approximation  to  the 
buildup  factor  is  given  by 


B(r)  =  1  +  a(r/X)k  , 

where  \  is  the  attenuation  length  and  a  may  be  shown  by  energy  conser¬ 
vation  to  be 

a  =  ys/[uar(k+1)]  , 

where 

Ug  =  total  scattering  coefficient  (ys  ”  og  +  1.022opp/E), 

Og  =  Compton  scattering  coefficient, 

opp  =  pair  production  attenuation  coefficient, 

E  =  gamma  ray  energy  in  MeV, 

Ua  =  total  absorption  coefficient  for  Compton  absorption, 
photoabsorption,  and  pair  production 
[pa  =  va  +  °ph  +  ®pp^2  *  1.022)/e]  , 

oa  =  Compton  absorption  coefficient, 

Oph  =  photoabsorption  attenuation  coefficient, 

T(k+1)  =  gamma  function  of  k+1, 

k  *  function  of  (source)  gamma  ray  energy. 

Parameter  k  is  of  order  unity  and  is  the  only  free  quantity  in  the 
expression  for  B. 

The  expression  for  B  was  fitted,  by  adjustment  of  k,  to  the  1-d 
detector  results  obtained  by  the  SLEDGE  code.  Twenty-five  detectors 
were  used,  each  one  attenuation  length  thick,  ltie  least-squares  fit  was 
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obtained  by  minimizing  the  sum  of  the  squares  of  the  errors  over  all 
detectors,  with  each  square  weighted  by  the  reciprocal  of  the  square  of 
the  fractional  deviation  for  the  detector.  Buildup  factor  parameters 
a  and  k  so  obtained  are  given  in  table  3  for  each  gamma  ray  energy. 
The  buildup  factor  predicted  by  these  values  of  a  and  k  has  a  probable 
error  of  not  more  than  1  percent  due  to  the  uncertainty  in  the  Monte 
Carlo  data.  Although  the  fractional  deviation  for  the  detectors  varies 
from  about  3  to  about  10  percent,  the  combined  fractional  deviation  (of 
all  25  detectors  fitted  at  once)  is  less  than  1  percent.  However,  this 
1-percent  error  estimate  does  not  include  the  possibility  of  systematic 
error  due  to  the  true  buildup 


factor  not  being  perfectly 

TABLE  3.  ENERGY  DEPOSITION 

BUILDUP  FACTOR 

represented  by  the  expression 

1  +  a(r/X)  .  This  error  esti¬ 

mate  does  describe  how  closely 
the  a  and  k  parameter  values. 

PARAMETERS 

Source  gamma  ray 

FOR  B  = 

=  1  +  a(r/X)* 

thus  obtained  by  fitting  Monte 
Carlo  data,  approximate  the  best 
fit  of  the  expression  to  the  true 

energy 

(MeV) 

a 

k 

buildup  factor. 

0.5 

1.534 

1.418 

As  to  possible  systemat¬ 

1.0 

1.162 

1.188 

ic  error,  it  appears  that  the  fit 
may  be  low  by  about  an  average  2 

1.5 

0.960 

1.  143 

percent  to  10  attenuation  lengths 

2.5 

0.800 

0.976 

range,  but  this  level  remains  un¬ 
certain  because  of  the  size  of 

3.5 

0.695 

0.903 

statistical  error  of  the  Monte 
Carlo  data.  Statistical  error  at 

5.0 

0.580 

0.880 

farther  ranges  is  even  greater 
and  permits  no  assessment  of  sys- 

7.0 

0  •  473 

0.891 

tematic  error. 


4.2  Ground  Effect  on  Buildup  Factor 

Analysis  of  the  2-d  detector  results  of  the  two  SLEDGE  runs, 
for  0.5-  and  1.5-MeV  gamma  ray  point  sources  at  200  m  aboveground, 
showed  a  distinct  depletion  of  energy  deposition  near  the  ground.  An 
exception  is  that  early  energy  deposition  is  enhanced  beneath  the  burst 
due  to  immediate  backscatter  from  the  ground.  However,  this  effect  will 
be  neglected  in  the  following  development  since  EMP  in  this  region  is 
not  important  to  typical  Army  electronic  systems  because  of  overwhelming 
blast  and  thermal  damage. 


For  observers  near  the  ground  and  for  a  range  such  that  the 
line  of  sight  to  the  burst  is  shallow,  the  buildup  factor  was  limited  in 
its  growth.  This  limitation  was  interpreted  that  the  ground  (although 


flat)  introduced  an  effective  horizon  because  the  ground  intercepts 
multiply  scattered  gamma  rays  that  would  otherwise  contribute.  The 
effective  horizon,  H,  is  approximately 


1/2 

H  =  3.0r/(9.0  +  r2/A2) 

for  horizontal  gamma  ray  incidence  near  the  ground. 

The  effective  horizon  for  horizontal  incidence,  h,  at  height  z 
aboveground  is  approximately 

h  =  Hez/X  . 


For  incidence  at  angle  0  to  the  horizontal,  the  effective  horizon  at  the 
observer  will  be  approximated  by  the  averaging  integral 


or 


h 


ez(r)/A 


[r( sin0)+Z] /A 


dr 


-  Aez/X  L H(sin0)/A 
sin  9  Le  \ 


where  z(r)  is  the  height  at  range  r  on  the  line  of  sight  from  the  burst 
to  the  observer  and  Z  is  the  observer  height  aboveground. 


The  ground- influenced  buildup  factor,  B(R,Z),  a  function  of 
observer  range  and  height,  is  asymptotically  limited  to  B(h),  that 
attained  within  effective  horizon  h.  Based  on  the  SLEDGE  2-d  results,  a 
good  approximation  for  the  ground- inf luenced  buildup  factor  is 

B(R,Z)  =  ( 1/2) [B(r)  +  X]  , 


where  B(r)  is  the  free-f ield  buildup  factor  for  range  r  from  the  source 
and 


X  =  1  +  1/(l/[B(r)  -  l]2  +  1/ {4  +  [B( h)  -  1] 2} ) 1/2  . 

4.3  Time  Dependence  of  Energy  Deposition  Rate 

The  time-dependent  energy  deposition  rate,  S,  at  distance  r, 
height  z,  and  time  t  contributed  by  multiply  scattered  gamma  rays  (that 
is,  excluding  uncollided  gamma  rays)  from  a  delta-function  (in  time) 
point  gamma  source  of  one  gamma  of  energy,  E,  is  written  as 

S(r,z,t)  *  U(r)N(r,z)T(r,z,t)  , 
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where 


U(r)  =  energy  deposition  from  uncollided  gamma  ray  number 

flux  (that  is,  first-scatter  deposition), 

N(r,z)  =  normalization  factor  (including  buildup  factor), 

T(r,z,t)  =  explicit  time  dependence. 

U(r)  is  given  by 

U(r)  =  Ee~r/V  (4irr2Xabs)  , 

where  E  is  the  gamma  ray  energy,  X  is  the  attenuation  length,  and  Xaks 
is  the  absorption  length.  Based  on  the  1.5-MeV  gamma  ray  source  SLEDGE 
results,  a  good  approximation  for  T  is 

T(r,z,t)  =  (t  +  tQ)m  e_pt  ' 

where 

tQ  =  108  s, 

m  =  — 1/ [ 1  +  f(r,z)r/X] , 

p  =  1.4  x  106  s'1. 

Function  f(r,z)  is  described  later. 

We  require  S  to  be  independent  of  z  at  t  =  0  since  the  ground 
affects  only  multiply  scattered  gamma  rays  (for  an  observer  above¬ 
ground)  ,  and  these  gamma  rays  always  arrive  after  t  =  0.  We  require 
also  the  following  normalization: 

J~Q  S(r,z,t)  dt  =  U(r)[S(r,z)  -  l]  , 


A 

where  B  is  the  ground- inf luenced  buildup  factor.  For  z  =  «®,  B  becomes 
free-field  buildup  factor  B(r).  From  this  normalization,  we  obtain 


N(r,z) 


T(r ,z,t) 


dt  =  B(r,z) 


1  . 


Let  I  be  defined  as  the  integral 


1  =  T(r,z,t)  dt 

-/.-(*  +  *o)’  e_Pt  «  • 
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By  straightforward  means#  the  integral  is  found  to  be 


I 


pto  . 
e  p 


cr(e)[l  -  p(e,pto)] 


where 

e  =  1  +  m, 

T(e)  =  gamma  function  of  e, 

P(e  pto)  =  incomPlete  gamma  function. ^ 

Hence,  normalization  factor  N(r,z)  is 

N(r,z)  =  [B(r,z)  -  l]/l  . 

Since  function  f  (r,z)  has  not  been  specified  yet,  variable  e, 
e  =  1  -  1/[1  +  f(r,z)r/X]  , 

remains  unconstrained.  We  use  it  to  satisfy  the  requirement  that 
S(r,z,t)  be  independent  of  z  at  t  =  0.  At  earliest  times,  the  transport 
is  essentially  free  field,  so  free-field  data  far  from  the  ground  should 

describe  the  t  =  0  transport.  By  examination  of  SLEDGE  data  for  detec¬ 

tors  far  from  the  ground  (virtually  free  field),  it  was  determined  that 
a  value  f  =  0.2  would  cause  S  to  fit  these  data  reasonably  well.  This 
value  is  called  fm  (that  is,  f„  =  0.2).  Let  the  corresponding  value  for 
e  be 

E„  -  1  “  V(l  +  f„rA)  • 

We  require  the  transport  near  the  ground  at  t  =  0,  to  be  the  same  as  the 
transport  far  (z  -*•  *)  from  the  ground,  also  at  t  =  0.  Thus, 

S(r,z=»,t=Q)  =  S(r,z,t=0)  , 


or 

[B( r)  -  13T(r,z=o»,t*0) 
I(r,z=») 

This  can  be  expanded  to 


[§( r ,z)  -  l]T(r,z,t=0) 
I(r,z) 


E«”l 

[B(r)  -  1] tp _ 

ePt°P  e"r(E»)[1  "  PCe—Pto)] 


pt-0 

e 


[§( r,z)  -  1  ]t 0  1 _ 

p“er(e)[l  -  P(e,pto) ] 
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5 Milton  Abramowitz  and  Irene  Stegrun,  Handbook  of  Mathematical 
Functions,  Dover  Publications,  Inc,,  New  York  (1965),  260, 


where  e  is  a  function  of  z.  If  we  define  G(e)  by 


G(e)  =  eept°p"er(e) [ 1  -  p(e,pt0)]  , 

then  this  may  be  written 

[B<r)  -  1]/G(eoo)  =  [£(r,z)  -  1]/G(e)j  , 
or 

G(e)  =  G(e»)[B(r,z)  -  l]/[B(r)  -  1]  . 

Thus,  given  r  and  z,  G(e)  may  be  determined  from  G (e^),  B(jr) ,  and 
B(r,z).  G(e)  may  then  be  inverted  to  yield  e.  Having  e  and  m  =  e  -  1, 
one  may  obtain 

I  =  G(e)to  , 

N(r,z)  =  [6(r ,z)  -  l]/I  , 

T(r ,z,t)  =  (t  +  to)me“pt  ' 

S(r ,z,t)  =  U(r)N(r,z)T(r,z,t)  . 


Inverting  G(e)  for  e  involves  tabulating  G  for  a  large  number  of  values 
of  e  and  then  using  a  table-lookup  algorithm  to  get  e  from  G. 

By  using  the  approach  just  described,  the  energy  deposition 
rate  may  be  obtained  for  any  r,z,t  in  such  a  way  that  it  evolves  in  time 
from  an  initial  free-field  distribution  to  a  progressively  more  ground- 
influenced  distribution.  Time  integrals  of  the  energy  deposition  rate 
are-  correctly  normalized  to  expected  ground- influenced  buildup  factors. 

These  results  are  scaled  to  first  order  to  source  gamma  ray 
energies  other  than  1.5  MeV  by  geometrically  scaling  r,  z,  and  t  by  the 
attenuation  length.  This  is  already  built  into  buildup  factors  B  and 
B  and  variable  m  ( and  e ) .  Only  the  time  remains  to  be  scaled.  If  sub¬ 
script  "1.5"  denotes  1.5-MeV  source  gamma  rays  and  if 

a  =  X,  /X (E)  , 

where  X(E)  is  the  attenuation  length  for  gamma  ray  energy  E,  then 
G(e)  =  (ato)1”^1*1  °p-cr<e>[l  -  p(c, potto)]  , 

I  *  G(e)(ato)m/a  , 


T(r,z,t)  »  (at  +  ato)me”pat  . 
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Variable  m  is  evaluated  for  A  =  A(E)  as 

m  =  -  1/Cl  +  fr/A(E)]  . 

The  uncollided  energy  deposition,  U,  becomes 

U(r)  =  Ee”r/x<E)/[4nr2Aaba(E) ]  . 

5.  RELATION  OF  COMPTON  CURRENT  TO  ENERGY  DEPOSITION  RATE 

5.1  Compton  Current  from  First  Collisions  of  Gamma  Rays 

For  point  source  gamma  rays  experiencing  their  first  collision 
in  air  having  density  1.11  kg/m3,  the  radial  Compton  current  produced  is 
simply  related  to  the  corresponding  energy  deposition: 

Jr/S  =  k0(E)  . 

Hie  theta  Compton  current,  Jg,  is  identically  zero  for  first  scatters. 
Factor  kg(E)  is  given  in  table  4  for  the  seven  gamma  ray  energies  con¬ 
sidered.  The  values  are  based  on  estimates6  of  the  mean  forward  range 
of  Compton  electrons,  including  multiple  scattering  of  the  Compton 
electrons  as  well.  Factor  kg(E)  equals  AaRc/(pAcE),  where  Aa  and  Ac 
are  the  Compton  absorption  and  collision  lengths,  p  is  the  air  density, 
R_  is  the  mean  forward  range  of  the  Compton  electrons,  and  E  is  the 
gamma  ray  energy. 

TABLE  4.  RATIO  OF  RADIAL  COMPTON  ELECTRON 
CURRENT  TO  ENERGY  DEPOSITION  FOR 
FIRST-COLLISION  GAMMA  RAYS  FROM 
A  POINT  SOURCE 


Source  gamma  ray 
energy 
(MeV) 

Ratio  k  t  E ) 
(electron-tn/MeV) 

0.5 

1.061 

1.0 

1.606 

1.5 

1.952 

2.5 

2.338 

3.5 

2.508 

5.0 

2.561 

7.0 

2.454 

6 William  T.  Hyatt,  Jr.,  Transmission  Factor  Effects  on  the  Average 
Forward  Range  of  Compton  Electrons,  Harry  Diamond  Laboratories  HDL-TH- 
80-10  (January  1980). 
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5.2  Radial  Compton  Current  from  Multiply  Scattered  Gamma  Rays 


For  multiply  scattered  gamma  rays,  the  radial  Compton  current 
produced  is  related  to  the  corresponding  energy  deposition  by 

Jj/S  =  k(r,t,E)  , 

where  r  is  the  distance  from  the  point  source  of  gamma  rays  and  t  is 
retarded  time.  The  SLEDGE  1-d  results  for  each  of  the  seven  gamma  ray 
energies  can  be  fitted  by 

k(r ,t,E)  *  ko (E) (o.91  +  0.03T.^/2)s~T  , 


where 
T  = 


1.0  x 


10“8pi(l.0  x 
[B(R)  -  1] 


108t)P2 

P4 


and 

R  =  r  +  X(p3  +  p5t/[r/A]P6)  . 

Factor  ko(E)  is  defined  in  section  5.1.  Buildup  factor  B  is  evaluated 
for  modified  range  R.  Optimum  fit  values  for  pj,  P2,  P3,  P4,  P5,  and 
p6  are  listed  in  table  5  for  the  seven  gamma  ray  energies. 


TABLE  5.  FIT  PARAMETERS  FOR  RATIO  OF  RADIAL  COMPTON  ELECTRON  CURRENT 
TO  ENERGY  DEPOSITION  RATE  FOR  MULTIPLY  SCATTERED  GAMMA  RAYS 


Source  gamma 
energy 
(MeV) 

ray 

P1 

P2 

P3 

P<, 

P5 

P6 

0.5 

3.880(7>a 

1.384 

7.5991-2) 

0.855 

2.793(7) 

-0.0847 

1.0 

1.783(7) 

1.040 

1. 304( -1) 

0.685 

0.817(7) 

+0.0524 

1.5 

1.  162(7) 

1.025 

7.  358 (-2) 

0.596 

0.537(7) 

-0.1179 

2.5 

0.964(7) 

0.974 

1.7651-2) 

0.679 

0.339(7) 

+0.0605 

3.5 

0.769(7) 

0.974 

1.679( -2 ) 

0.758 

0.247(7) 

+0.1410 

5.0 

0.358(7) 

1.230 

7.253(-2) 

0.947 

0.523(7) 

-0.1091 

7.0 

0.835(6) 

1.865 

7. 869 (-2) 

1.561 

1.336(7) 

-0.0886 

Quantity  k(r,t,E)  decays  more  slowly  with  increasing  range. 
Also,  in  theory  one  would  expect  k  at  time  t  5  0  to  be  exactly  ko(E) 
since  the  earliest  scattered  gamma  rays  to  arrive  at  any  point  are  those 
that  have  experienced  a  single  grazing  scatter  and  thus  have  nearly  the 
same  energy,  E.  However,  the  Monte  Carlo  studies  show  that 

k(r ,t,E)  *  kO(E)(o.91  +  0.03El/z) 

when  averaged  over  the  first  time  bin,  which  is  5  ns  wide.  This  resolu¬ 
tion  of  the  time  history  appears  to  be  adequate  since  most  gamma  ray 
output  from  a  nuclear  burst  is  emitted  over  about  10  ns. 

5.3  Theta  Compton  Current  from  Multiply  Scattered  Gamma  Rays 

Hie  theta  Compton  current  from  a  point  gamma  ray  source  near 
the  earth  is  secondary  to  the  radial  Compton  current  in  magnitude  and  in 
efficiency  as  a  generator  of  EMP.  In  the  SLEDGE  2-d  results,  variance 
was  larger  for  the  theta  current  component — often  so  large  that  the  sign 
of  the  component  was  doubtful.  Nevertheless,  some  significant  trends 
were  clearly  observed  at  late  times  near  the  ground.  Hie  ratio  of  theta 
Compton  current  to  energy  deposition  rate  (for  an  air  density  of  1.11 
kg/m3 )  is  approximately 


Je/S  =  1.5  x  105te'0,02z  , 

where  t  is  retarded  time  in  seconds  and  z  is  height  aboveground  in 
meters.  Hius,  the  theta  component  is  strongest  near  the  ground  and 
decays  (about)  exponentially  with  increasing  height. 

This  expression  for  Jg/S  is  not  strongly  dependent  upon  source 
gamma  ray  energy  to  a  first  approximation  since  the  late-time  currents 
derive  from  low-energy  multiply  scattered  gamma  rays. 


6.  COMPARISON  WITH  PREVIOUS  RESULTS 

6.1  Buildup  Factor  for  Energy  Deposition 

The  LEMP  code  prompt  gamma  ray  drivers  (H.  J.  Langley,  Mission 
Research  Corp.)  include  buildup  factor  curve  fits  for  energy  deposition 
and  radial  Compton  current  for  three  source  gamma  ray  energies  (0.5, 
1.5,  and  5.0  MeV).  No  account  is  taken  of  the  influence  of  the  ground. 
The  energy  deposition  buildup  factor  obtained  for  air  was  based  on 
published  data  for  water. 

Recently,  Malik,  Cashwell,  and  Schrandt7  at  the  Los  Alamos 
Scientific  Laboratory  ( LASL)  obtained  curve  fits  to  buildup  factors  and 

7J.  5.  Malik,  E.  D.  Cashwell,  and  R.  G.  Schrandt,  The  Time  Dependence 
of  the  Compton  Current  and  Energy  Deposition  from  Scattered  Gamma  Rays, 
Los  Alamos  Scientific  Laboratory,  EM,  LA-7386-MS  (July  1978), 
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time  dependence  of  energy  deposition  and  Compton  current  from  scattered 
gamma  rays  in  homogeneous  isotropic  air.  These  curve  fits  to  their 
Monte  Carlo  data  were  obtained  for  gamma  ray  energies  of  1,  1.5,  2,  3, 
4,  and  6  MeV  to  a  range  of  1600  m.  Air  density  was  1.1  mg/cm3. 

The  LEMP  and  LASL  buildup  factors  are  compared  with  HEMP  free- 
f ield  and  ground-influenced  buildup  factors  ( for  energy  deposition)  in 
figures  1  to  6.  The  NEMP  ground- inf luenced  buildup  factor  is  for  a 
burst  and  an  observer  both  at  a  0-m  height  aboveground.  Figures  1  to  3 
show  results  as  a  function  of  range  in  units  of  attenuation  lengths, 
from  0  to  25  attenuation  lengths.  These  results  exclude  the  effect  of 
any  disagreement  in  attenuation  length  among  the  LEMP,  LASL,  and  NEMP 
drivers  and  compare  directly  the  functional  forms  used.  Figures  1  to  3 
show,  respectively,  0.5-,  1.5- ,  and  5.0-MeV  source  gamma  ray  buildup 
factors.  The  0.5-MeV  energy  is  below  the  1.0-MeV  lower  limit  for  the 
LASL  transport  data  and  shows  that  the  LASL  drivers  are  in  substantial 
error  below  the  designated  range.  The  LEMP  and  NEMP  curves  are  in  good 
agreement  to  about  four  attenuation  lengths  range,  after  which  the  LEMP 
curve  is  higher.  For  0.5  MeV,  the  disagreement  becomes  very  large.  Hie 
quality  of  agreement  between  LEMP  and  NEMP  at  closer  ranges,  however, 
tends  to  confirm  their  validity  at  such  ranges.  The  LASL  curves,  on  the 
other  hand,  seem  too  low  except  for  the  5.0-MeV  results  (fig.  3).  The 
low-energy  disagreement  among  LEMP,  NEMP,  and  LASL  is  ameliorated  by  the 
fact  that,  for  a  realistic  gamma  ray  energy  spectrum,  higher  energy 
gamma  rays  predominate  at  substantial  ranges  because  of  their  greater 
penetrating  ability. 


Figure  1.  Energy  deposition  buildup  factors 
compared  for  0.5-MeV  gamma  rays,  for  range 
in  attenuation  lengths. 
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Figure  2.  Energy  deposition  buildup  factors 
compared  for  1.5-MeV  gamma  rays,  for  range 
in  attenuation  lengths. 
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Figure  3.  Energy  deposition  buildup  factors 
compared  for  5.0-MeV  gamma  rays,  for  range 
in  attenuation  lengths. 
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Figure  6.  Energy  deposition  buildup  factors 
compared  for  5.0-MeV  gamma  rays,  for  range 
in  meters. 


Figures  4  to  6  include  the  attenuation  length,  but  still  appear 
similar  to  figures  1  to  3.  Careful  comparison  of  SLEDGE  transport 
results  with  the  NEMP  curve  fit  shows  that  the  NEMP  curve  fit  is  about 
10  percent  higher  over  the  first  attenuation  length.  But  the  NEMP  curve 
fit  is  found  to  be  accurate  to  a  few  percent  at  greater  ranges. 

The  NEMP  ground- inf luenced  buildup  factor  in  figures  1  to  6  is 
as  much  a  factor  of  two  smaller  than  the  corresponding  NDiP  free-field 
buildup  factor.  The  difference  increases  monotonically  with  range,  more 
slowly  for  more  energetic  gamma  rays,  to  a  difference  of  as  much  as  a 
factor  of  two. 

6.2  Time  Dependence  of  Btergy  Deposition  Rate 

The  time  dependence  of  the  free-field  NEMP,  LEMP,  and  LASL 
energy  deposition  rate  drivers  is  shown  in  figures  7  to  12  for  three 
source  gamma  ray  energies  (0.5,  1.5,  and  5.0  MeV)  and  two  ranges  from 
the  burst  (300  and  3000  m) .  The  three  curves  are  normalized  to  a  time 
integral  of  unity,  removing  any  disagreement  in  attenuation  length  and 
buildup  factor.  A  fourth  curve  represents  the  ground- influenced  NEMP 
energy  deposition  rate  for  burst  and  observer  on  the  ground.  This  curve 
is  normalized  to  the  free-field  case,  so  the  difference  between  this 
curve  and  the  free-field  NB1P  curve  is  a  diminution  of  the  energy  depo¬ 
sition  rate  due  to  the  presence  of  the  ground. 
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Figure  7.  Time  dependence  of  energy  deposition 
buildup  factors  compared  at  300-m  range  for 
0.5-MeV  gamma  rays. 


Figure  8.  Time  dependence  of  energy  deposition 
buildup  factors  compared  at  300-m  range  for 
1.5-MeV  gamma  rays. 
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Figure  11.  Time  dependence  of  energy  deposition 
buildup  factors  compared  at  3000-m  range  for 
1.5-MeV  gamma  rays. 


Figure  12 .  Time  dependence  of  energy  deposition 
buildup  factors  compared  at  3000-m  range  for 
5.0-MeV  gamma  rays. 


Figures  7  to  9  show  the  energy  deposition  rate  at  a  300-m  range 
for  0.5-,  1.5-,  and  5.0-MeV  source  gamma  ray  energies.  The  agreement 
between  the  NEMP  and  LEMP  drivers  is  good.  Since  the  LEMP  driver  was 
based  on  close-in  data,  we  expect  it  to  be  accurate  at  a  300-ra  range. 
The  extremely  close  agreement  at  1-ns  retarded  time  is  particularly 
encouraging  and  tends  to  confirm  the  correctness  of  the  first-scatter 
part  of  the  NEMP  driver  •  The  LASL  curve  seems  to  be  concentrated  too 
much  toward  early  times,  except  for  the  0.5-MeV  curve. 

Figures  10  to  12  show  the  energy  deposition  rate  at  a  3000-m 
range  for  0.5-,  1.5-,  and  5.0-MeV  source  gamma  ray  energies.  At  this 
deep  penetration,  we  expect  the  LEMP  curves  to  be  concentrated  too  much 
toward  early  times  in  comparison  with  the  NEMP  curve — and  they  are.  The 
NEMP  and  LASL  curves  agree  well  only  for  1.5  MeV  (fig.  11).  For  5.0  MeV 
(fig.  12),  the  LASL  curve  is  too  close  to  the  LEMP  curve.  Data  shown  in 
the  LASL  report7  confirm  a  tendency  of  that  curve  fit  to  overestimate 
the  energy  deposition  rate  earlier  than  a  fraction  of  a  microsecond,  for 
deep  penetrations  of  about  10  attenuation  lengths. 

In  all  figures  shown,  the  NEMP  ground- influenced  energy  depo¬ 
sition  rate  begins  equal  to  the  NEMP  free-field  curve  and  later  falls 
below  by  as  much  as  a  factor  of  about  three  for  the  3000-m  curve.  For 
the  300-m  curve,  little  difference  is  noted. 

6.3  Time  Dependence  of  Radial  Compton  Current 

The  ratio  of  radial  Compton  electron  current  to  energy  depo¬ 
sition  rate  is  shown  as  a  function  of  time  for  the  NEMP,  LEMP,  and  LASL 
drivers  in  figures  13  to  18.  Figures  13  to  15  are  for  gamma  ray  ener¬ 
gies  of  0.5,  1.5,  and  5.0  MeV  at  a  range  of  300  m  from  the  burst. 
Figures  16  to  18  are  for  the  same  gamma  ray  energies  at  a  range  of  3000 
m.  The  units  of  the  ratios  plotted  are  in  Compton  electron-meter  per 
MeV.  A  factor  of  0.8  has  been  separately  applied  to  the  LEMP  and  LASL 
drivers  to  account  for  multiple  scattering  of  Compton  electrons.  This 
effect  is  already  included  in  the  NEMP  driver.  As  for  energy  deposition 
rate,  we  expect  the  LEMP  data  to  be  more  accurate  at  300  m  than  at  3000 
m.  Good  agreement  is  observed  among  all  three  curves  (NEMP,  LEMP,  and 
LASL)  at  300  m  for  each  gamma  ray  energy,  with  a  few  exceptions.  The 
NEMP  curve  is  probably  too  high  after  about  1  (is,  although  the  signifi¬ 
cance  of  the  discrepancy  is  slight  in  terms  of  the  EMP  that  would  be 
produced.  Also,  the  LASL  early  time  value  for  0.5  MeV  is  high  by  more 
than  a  factor  of  two. 


7J.  S .  Malik,  E.  D.  Castnrell,  and  R.  G.  Schrandt,  The  Time  Dependence 
of  the  Compton  Current  and  Energy  Deposition  from  Scattered  Gamma  Rays , 
Los  Alamos  Scientific  Laboratory,  MM,  LA-7 386-MS  (July  1978). 
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Figure  14.  Time  dependence  compared  for  ratio 
of  radial  Compton  electron  current  to  energy 
deposition  rate,  at  300-m  range  for  1.5-MeV 
gamma  rays. 
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Agreement  is  worse  at  3000  m  for  each  gamma  ray  energy.  ttie 
NEMP  and  LEMP  curves  are  in  fair  agreement  at  early  times,  whereas  the 
LASL  curve  is  too  high  except  for  the  5.0-MeV  curve.  The  LEMP  curve 
decays  much  too  fast,  as  expected  at  deep  penetration  where  multiply 
scattered  gamma  rays  retain  more  of  their  collimation  than  at  shallow 
penetration.  (The  LEMP  driver  is  based  on  shallow  penetration  data). 
The  LASL  curve,  on  the  other  hand,  decays  too  slowly.  Comparison  of  the 
NEMP  curves  with  the  Monte  Carlo  data  that  they  were  intended  to  fit 
reveals  that  the  curves  are  usually  somewhat  too  high  after  about  1  us, 
up  to  a  factor  of  about  two.  However,  considering  the  steepness  of  the 
curve,  this  error  is  not  serious.  No  doubt  a  more  complicated 
functional  form  would  approximate  the  Monte  Carlo  data  more  faithfully. 

Fran  these  comparisons,  LEMP  drivers  would  underestimate  the 
Compton  current  at  3000  m  by  about  a  factor  of  four  at  100  ns  after  the 
gamma-dot  peak,  for  any  of  the  gamma  ray  energies  considered.  The  LASL 
drivers  should  not  be  applied  to  gamma  ray  energies  less  than  1.0  MeV, 
and  they  also  significantly  overestimate  the  Compton  current  at  3000  m 
for  1.5  MeV. 

No  ground  effect  on  the  ratio  of  Compton  current  to  energy 
deposition  rate  was  identified. 


7.  CONCLUDING  REMARKS 

This  report  presents  improved  calculations  of  EMP  drivers  due  to  a 
gamma  ray  point  source  in  air  over  ground  and  analytic  approximations  to 
these  EMP  drivers.  These  were  done  for  seven  source  gamma  ray  energies 
from  0.5  to  7.0  MeV.  The  EMP  driver  calculations  were  based  on  Monte 
Carlo  simulation  of  the  gamma  ray  transport  using  the  SLEDGE  computer 
code  developed  for  this  study.  Notable  features  of  the  simulation 
included  tracking  to  25  attenuation  lengths  (mean  free  paths)  penetra¬ 
tion,  presence  of  the  earth,  use  of  correlated  sampling  to  estimate 
Compton  electron  currents  more  accurately,  and  extensive  use  of  biasing 
and  importance  functions.  The  theoretical  basis  of  the  SLEDGE  code  is 
discussed.  Significant  aspects  of  the  Monte  Carlo  results  also  are 
discussed. 

The  analytic  approximations  to  the  EMP  drivers  are  described  in 
detail,  including  (1)  the  energy  deposition  buildup  factor  (both  free 
field  and  ground  influenced),  (2)  time  dependence  of  energy  deposition 
as  a  function  of  both  range  from  the  burst  and  nearness  to  the  ground, 
(3)  ratio  of  radial  Compton  electron  current  to  energy  deposition  rate, 
and  (4)  theta  Compton  electron  current  due  to  the  presence  of  the 
ground.  These  approx imat' 'ns  have  been  incorporated  into  an  EMP  driver 
software  package  used  with  the  NEMP  near-surface-burst  EMP  environment 
prediction  code.  Comparisons  are  made  with  two  previous  EMP  driver 


packages,  one  used  in  the  LEMP  surface-burst  EMP  environment  prediction 
code  and  the  other  originated  by  LASL  scientists  for  high-altitude  burst 
problems*  Free- field  agreement  with  the  LEMP  drivers  was  generally  good 
at  short  ranges.  Free-field  agreement  with  the  LASL  drivers  was  better 
at  higher  gamma  ray  energies.  The  new  results  predict  ground  effects  as 
well,  which  were  not  considered  by  the  LEMP  and  LASL  researchers. 
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ATTN  TRAVIS  L.  SIMPSON 

LUBBOCK,  TX  79417 

UNITED  TECHNOLOGIES  CORP 
HAMILTON  STANDARD  OI VIS ION 
BRADLEY  INTERNATIONAL  AIRPORT 
ATTN  CHIEF  ELEC  DESIGN 

WINDSOR  LOCKS,  CT  06069 

WEST IN GHO USE  ELECTRIC  CORPORATION 
ADVANCED  ENERGY  SYSTEMS  DIV 
P.O.  BOX  10864 
ATTN  TECH  LIB 
PITTSBURGH,  PA  15236 

US  ARMY  ELECTRONICS  RESEARCH 
6  DEVELOPMENT  COMMAND 
ATTN  TECHNICAL  DIRECTOR,  DRDEL-CT 

HARRY  DIAMOND  LABORATORIES 

ATTN  CO/  TD/TSO/DI  VIS  ION  DIRECTORS 

ATTN  RECORD  COPY,  81200 

ATTN  KDL  LIBRARY  81100  (3  COPIES) 

ATTN  HOL  LIBRARY  (WOODBRIDGB)  81100 

ATTN  TECHNICAL  REPORTS  BRANCH,  81300 

ATTN  CHAIRMAN,  EDITORIAL  COMMITTEE 
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CHIEF, 

21000 

ATTN 

CHIEF, 

22000 

ATTN 

CHIEF, 

22100 

(3 

COPIES) 

ATTN 

CHIEF, 

22300 

ATTN 

CHIEF, 

22800 

ATTN 

CHIEF, 

22900 

ATTN 

CHIEF, 

13300 

ATTN 

CHIEF, 

21100 

(3 

COPIES) 

ATTN 

CHIEF, 

21200 

ATTN 

CHIEF, 

21300 

(5 

COPIES ) 

ATTN 

CHIEF, 

21400 

(2 

COPIES) 

ATTN 

CHIEF, 

21500 

ATTN  BALICKI,  F.  W.,  20240 

ATTN  WIMENITZ,  P.  N.,  20240 

ATTN  TALLERICO,  A.,  47400 

ATTN  BIXBY,  R.,  22900 

ATTN  WYATT,  W.  T. ,  21300  (20  COPIES) 
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